# Replication package for 
# "The Economic Leverage of International Organizations in Interstate Disputes"
# Johannes Karreth
# June 30, 2017
# jkarreth@ursinus.edu

# This file: 24_crises_coll_analyses.R
# Purpose: Analyses of major clashes or war during crises, with contemporaneous crises collapsed

# Note: must have read in data before, using 21_crises_readdata.R

# File structure:
# 1. Crises collapsed to dyads with largest powers
# 2. Crises collapsed to crises with (most) violent outcomes

# setwd("...")

library("rstan")
library("rstanarm")
library("texreg")
library("ggplot2")
library("xtable")
library("gridExtra")
library("loo")
library("lme4")
library("dplyr")
library("reshape2")
library("BMA")

# Source in functions
source("Functions/MCMClogit.fd.mat.R")
source("Functions/theme_jk.R")
source("Functions/plotBMA.R")
source("Functions/MCMC_roc_prc.r")
source("Functions/geom_flat_violin.R")

# Unit of analysis: crisis (collapsed to dyads with largest powers)

# Model 1

m1_eq <- viol.majclash ~ igo_lev3_count_use + highgrav + territory + jointpol7 + rivalry_th + absidealdiff + atopally
m1_mf <- model.frame(m1_eq, data = dat.noclus)

m1 <- glm(m1_eq, data = m1_mf, family = binomial(link = "logit"))

m1_stanglm <- stan_glm(formula = m1_eq,
                       data = (m1_mf),
                       # family = binomial(link = "logit"),
                       # prior = student_t(df = 7, location = 0, scale = 2.5), 
                       # prior_intercept = normal(0, 10),
                       family = binomial(link = "logit"),
                       prior = normal(location = 0, scale = 10), 
                       prior_intercept = normal(0, 10),
                       seed = 20170207,
                       cores = 4)

saveRDS(m1_stanglm, file = "Output_MCMC/crises_noclus_m1_stanglm.RDS")

# Table

m1_tab <- m1_stanglm$stan_summary[, c(1, 3)]
m1_tab <- data.frame(m1_tab[c(2:8, 1, 10), ])
m1_tab$varname <- c("IGOs with high leverage", 
                    "Existential threat", 
                    "Territorial dispute", 
                    "Joint democracy", 
                    "Strategic rivalry",
                    "UNGA ideal point difference",
                    "Allies",
                    "Intercept",
                    "Log-posterior density")
names(m1_tab) <- c("Mean", "Standard deviation", "Variable")
m1_tab <- m1_tab[, c(3, 1, 2)]

m1_xtab <- xtable(m1_tab, digits = 2)
print(m1_xtab, file = "Output_Tables-and-Figures/crises_noclus_m1_table.tex", include.rownames = FALSE)

# First differences

m1_mm <- model.matrix(m1_eq, data = dat.noclus)

temp_sims <- as.matrix(m1_stanglm)
m1.fd <- MCMClogit.fd.mat(model_matrix = m1_mm, mcmc_out = temp_sims, 
                          credint = c(0.025, 0.975), 
                          percentiles = c(0.1, 0.9),
                          full_sims = TRUE)

m1.changes <- matrix(NA, ncol = 2, nrow = ncol(m1_mm))
for (i in 2:ncol(m1_mm)){
  m1.changes[i, ] <- ifelse(length(unique(m1_mm[, i])) == 2 & range(m1_mm[, i]) == c(0, 1), c(0, 1), 
                            quantile(m1_mm[, i], probs = c(0.1, 0.9)))
}

# > m1.changes
# [,1]  [,2]
# [1,]    NA    NA
# [2,] 0.000 5.000
# [3,] 0.000 1.000
# [4,] 0.000 1.000
# [5,] 0.000 1.000
# [6,] 0.000 1.000
# [7,] 0.143 4.055
# [8,] 0.000 1.000

fd.dat <- data.frame(m1.fd)
fd.dat$variable_label <- c("IGOs with high leverage\n(from 0 to 5)", 
                           "Existential threat\n(no to yes)", 
                           "Territorial dispute\n(no to yes)", 
                           "Joint democracy\n(no to yes)", 
                           "Strategic rivalry\n(no to yes)",
                           "UNGA ideal point difference\n(0.1 to 4)",
                           "Allies \n(no to yes)")

fd_long <- melt(m1.fd)
fd_long$Varpos <- ifelse(fd_long$Var2 == "igo_lev3_count_use", 1,
                         ifelse(fd_long$Var2 == "highgrav", 2, 
                                ifelse(fd_long$Var2 == "territory", 3,
                                       ifelse(fd_long$Var2 == "jointpol7", 4,
                                              ifelse(fd_long$Var2 == "rivalry_th", 5,
                                                     ifelse(fd_long$Var2 == "absidealdiff", 6, 
                                                            ifelse(fd_long$Var2 == "atopally", 7,
                                                                   NA)))))))

fd_long$variable_label <- ifelse(fd_long$Var2 == "igo_lev3_count_use", "IGOs with high leverage\n(from 0 to 5)",
                                 ifelse(fd_long$Var2 == "highgrav", "Existential threat\n(no to yes)", 
                                        ifelse(fd_long$Var2 == "territory", "Territorial dispute\n(no to yes)",
                                               ifelse(fd_long$Var2 == "jointpol7", "Joint democracy\n(no to yes)",
                                                      ifelse(fd_long$Var2 == "rivalry_th", "Strategic rivalry\n(no to yes)",
                                                             ifelse(fd_long$Var2 == "absidealdiff", "UNGA ideal point difference\n(0.1 to 4)", 
                                                                    ifelse(fd_long$Var2 == "atopally", "Allies \n(no to yes)",
                                                                           NA)))))))

# ROPE: use SE of the ratio of 1s

r_n <- length(m1_mf$viol.majclash)
r_p <- sum(m1_mf$viol.majclash) / r_n
r_se <- sqrt((r_p * (1 - r_p)) / r_n)
ROPE <- r_se

fd.p <- ggplot(data = fd_long, aes(x = reorder(variable_label, -Varpos), y = value)) + 
  scale_y_continuous(labels = function(x) x*100) +
  geom_flat_violin(color = NA, fill = "black") +
  coord_flip() + 
  theme_classic() + 
  ylab("Percentage point change in the probability of major clashes\nduring a crisis as each predictor changes as indicated") + 
  xlab("") + 
  theme(axis.ticks = element_blank(), axis.text = element_text(color = "black"), axis.line = element_blank()) + 
  annotate(geom = "text", y = -0.05, x = 0.15, label = "Major clashes less likely", hjust = 1, size = 3) +
  annotate(geom = "text", y = 0.05, x = 0.15, label = "Major clashes more likely", hjust = 0, size = 3) +
  annotate(geom = "segment", y = 0, yend = -0.7, x = 0, xend = 0, arrow = arrow(ends = "last", type = "closed"), size = 1) +
  annotate(geom = "segment", y = 0, yend = 0.7, x = 0, xend = 0, arrow = arrow(ends = "last", type = "closed"), size = 1) +
  annotate(geom = "rect", ymin = 0 - ROPE, ymax = 0 + ROPE, xmin = 0, xmax = Inf, fill = "gray", color = NA, alpha = 0.5) +
  geom_segment(y = 0, yend = 0, x = 0, xend = 10.5, color = "black", linetype = "solid", size = 0.33) + 
  annotate(geom = "text", y = -0.15, x = 5.25, hjust = 1, vjust = 0.5, label = "Region of practical\nequivalence (ROPE)\n[-2, 2]", color = "darkgray", size = 3) +
  annotate(geom = "segment", y = -0.15, yend = -1 * ROPE, x = 5.25, xend = 5.3, size = 0.25, color = "darkgray") +
  annotate(geom = "text", y = 0.1, x = 7.25, hjust = 0, vjust = 0.5, label = "% of posterior distribution\nto the left of ROPE", size = 3) +
  annotate(geom = "segment", y = 0.1, yend = -0.2, x = 7.25, xend = 7.15, size = 0.25) +
  annotate(geom = "text", y = 0.35, x = 4.35, hjust = 0.5, vjust = 0, label = "% of posterior distribution\nto the right of ROPE", size = 3) +
  annotate(geom = "segment", y = 0.3, yend = 0.2, x = 4.35, xend = 4.05, size = 0.25) +
  annotate(geom = "text", y = -0.45, x = 6.5, hjust = 0.5, vjust = 0, label = "Mean estimated change", size = 3) +
  annotate(geom = "segment", y = -0.5, yend = -0.3, x = 6.65, xend = 6.9, size = 0.25)

# % < 0

m1.fd_out0 <- apply(m1.fd, 2, function(x) ifelse(median(x) < 0, sum(x < 0) / length(x), sum(x > 0) / length(x)))
m1.fd_outROPE <- apply(m1.fd, 2, function(x) ifelse(median(x) < 0, sum(x < 0 - ROPE) / length(x), sum(x > 0 + ROPE) / length(x)))

fd_changelab <- ifelse(apply(m1.fd, 2, median) < 0, "Decrease by", "Increase by")
fd_changelab_short <- ifelse(apply(m1.fd, 2, median) < 0, "-", "+")
fd_changetext <- paste(fd_changelab, abs(round(apply(m1.fd, 2, median) * 100)), "points", sep = " ")
fd_changetext_short <- paste(fd_changelab_short, abs(round(apply(m1.fd, 2, median) * 100)), sep = "")

m1.fd_annotate <- data.frame(x = apply(m1.fd, 2, median), 
                             y = rev(unique(fd_long$Varpos)), 
                             out0 = paste(round(m1.fd_out0 * 100, digits = 1), "%", sep = ""), 
                             outROPE = paste(round(m1.fd_outROPE * 100, digits = 1), "%", sep = ""),
                             change = fd_changetext_short)
fd.p <- fd.p + geom_text(data = m1.fd_annotate, aes(y = x, x = y, label = outROPE), color = "white", nudge_x = 0.075, size = 3) + 
  geom_text(data = m1.fd_annotate, aes(y = x, x = y, label = change), color = "black", nudge_x = -0.1, size = 4)

ggsave(fd.p, file = "Output_Tables-and-Figures/crises_noclus_m1_fd.pdf", width = 6, height = 8)

# Unit of analysis: Crises collapsed to crises with (most) violent outcomes

# Model 1

m1_eq <- viol.majclash ~ igo_lev3_count_use + highgrav + territory + jointpol7 + rivalry_th + absidealdiff + atopally
m1_mf <- model.frame(m1_eq, data = dat_nomult)

m1 <- glm(m1_eq, data = m1_mf, family = binomial(link = "logit"))

m1_stanglm <- stan_glm(formula = m1_eq,
                       data = (m1_mf),
                       # family = binomial(link = "logit"),
                       # prior = student_t(df = 7, location = 0, scale = 2.5), 
                       # prior_intercept = normal(0, 10),
                       family = binomial(link = "logit"),
                       prior = normal(location = 0, scale = 10), 
                       prior_intercept = normal(0, 10),
                       seed = 20170207,
                       cores = 4)

saveRDS(m1_stanglm, file = "Output_MCMC/crises_nomult_hiviol_m1_stanglm.RDS")

# Table
m1_tab <- m1_stanglm$stan_summary[, c(1, 3)]
m1_tab <- data.frame(m1_tab[c(2:8, 1, 10), ])
m1_tab$varname <- c("IGOs with high leverage", 
                    "Existential threat", 
                    "Territorial dispute", 
                    "Joint democracy", 
                    "Strategic rivalry",
                    "UNGA ideal point difference",
                    "Allies",
                    "Intercept",
                    "Log-posterior density")
names(m1_tab) <- c("Mean", "Standard deviation", "Variable")
m1_tab <- m1_tab[, c(3, 1, 2)]

m1_xtab <- xtable(m1_tab, digits = 2)
print(m1_xtab, file = "Output_Tables-and-Figures/crises_nomult_hiviol_m1_table.tex", include.rownames = FALSE)

# First differences

m1_mm <- model.matrix(m1_eq, data = dat_nomult)

temp_sims <- as.matrix(m1_stanglm)
m1.fd <- MCMClogit.fd.mat(model_matrix = m1_mm, mcmc_out = temp_sims, 
                          credint = c(0.025, 0.975), 
                          percentiles = c(0.1, 0.9),
                          full_sims = TRUE)

m1.changes <- matrix(NA, ncol = 2, nrow = ncol(m1_mm))
for (i in 2:ncol(m1_mm)){
  m1.changes[i, ] <- ifelse(length(unique(m1_mm[, i])) == 2 & range(m1_mm[, i]) == c(0, 1), c(0, 1), 
                            quantile(m1_mm[, i], probs = c(0.1, 0.9)))
}

# > m1.changes
# [,1]  [,2]
# [1,]    NA    NA
# [2,] 0.000 5.000
# [3,] 0.000 1.000
# [4,] 0.000 1.000
# [5,] 0.000 1.000
# [6,] 0.000 1.000
# [7,] 0.143 4.055
# [8,] 0.000 1.000

fd.dat <- data.frame(m1.fd)

fd_long <- melt(m1.fd)
fd_long$Varpos <- ifelse(fd_long$Var2 == "igo_lev3_count_use", 1,
                         ifelse(fd_long$Var2 == "highgrav", 2, 
                                ifelse(fd_long$Var2 == "territory", 3,
                                       ifelse(fd_long$Var2 == "jointpol7", 4,
                                              ifelse(fd_long$Var2 == "rivalry_th", 5,
                                                     ifelse(fd_long$Var2 == "absidealdiff", 6, 
                                                            ifelse(fd_long$Var2 == "atopally", 7,
                                                                   NA)))))))

fd_long$variable_label <- ifelse(fd_long$Var2 == "igo_lev3_count_use", "IGOs with high leverage\n(from 0 to 5)",
                                 ifelse(fd_long$Var2 == "highgrav", "Existential threat\n(no to yes)", 
                                        ifelse(fd_long$Var2 == "territory", "Territorial dispute\n(no to yes)",
                                               ifelse(fd_long$Var2 == "jointpol7", "Joint democracy\n(no to yes)",
                                                      ifelse(fd_long$Var2 == "rivalry_th", "Strategic rivalry\n(no to yes)",
                                                             ifelse(fd_long$Var2 == "absidealdiff", "UNGA ideal point difference\n(0.1 to 4)", 
                                                                    ifelse(fd_long$Var2 == "atopally", "Allies \n(no to yes)",
                                                                           NA)))))))

# ROPE: use SE of the ratio of 1s

r_n <- length(m1_mf$viol.majclash)
r_p <- sum(m1_mf$viol.majclash) / r_n
r_se <- sqrt((r_p * (1 - r_p)) / r_n)
ROPE <- r_se

fd.p <- ggplot(data = fd_long, aes(x = reorder(variable_label, -Varpos), y = value)) + 
  scale_y_continuous(labels = function(x) x*100) +
  # scale_x_discrete(labels = rev(fd.dat$variable_label)) +
  # geom_rect(aes(ymin = 0 - ROPE, ymax = 0 + ROPE, xmin = 0, xmax = Inf), col = "gray") +
  # geom_segment(y = 0, yend = 0, x = 0, xend = 10.5, color = "black", linetype = "solid", size = 0.33) + 
  geom_flat_violin(color = NA, fill = "black") +
  coord_flip() + 
  theme_classic() + 
  ylab("Percentage point change in the probability of major clashes\nduring a crisis as each predictor changes as indicated") + 
  xlab("") + 
  # scale_y_discrete(labels = rev(fd.dat$variable_label), expand = c(0, 0.6)) +
  theme(axis.ticks = element_blank(), axis.text = element_text(color = "black"), axis.line = element_blank()) + 
  annotate(geom = "text", y = -0.05, x = 0.15, label = "Major clashes less likely", hjust = 1, size = 3) +
  annotate(geom = "text", y = 0.05, x = 0.15, label = "Major clashes more likely", hjust = 0, size = 3) +
  annotate(geom = "segment", y = 0, yend = -0.7, x = 0, xend = 0, arrow = arrow(ends = "last", type = "closed"), size = 1) +
  annotate(geom = "segment", y = 0, yend = 0.7, x = 0, xend = 0, arrow = arrow(ends = "last", type = "closed"), size = 1) +
  annotate(geom = "rect", ymin = 0 - ROPE, ymax = 0 + ROPE, xmin = 0, xmax = Inf, fill = "gray", color = NA, alpha = 0.5) +
  geom_segment(y = 0, yend = 0, x = 0, xend = 10.5, color = "black", linetype = "solid", size = 0.33) + 
  annotate(geom = "text", y = -0.15, x = 5.25, hjust = 1, vjust = 0.5, label = "Region of practical\nequivalence (ROPE)\n[-2, 2]", color = "darkgray", size = 3) +
  annotate(geom = "segment", y = -0.15, yend = -1 * ROPE, x = 5.25, xend = 5.3, size = 0.25, color = "darkgray") +
  annotate(geom = "text", y = 0.1, x = 7.25, hjust = 0, vjust = 0.5, label = "% of posterior distribution\nto the left of ROPE", size = 3) +
  annotate(geom = "segment", y = 0.1, yend = -0.2, x = 7.25, xend = 7.15, size = 0.25) +
  annotate(geom = "text", y = 0.35, x = 4.35, hjust = 0.5, vjust = 0, label = "% of posterior distribution\nto the right of ROPE", size = 3) +
  annotate(geom = "segment", y = 0.3, yend = 0.2, x = 4.35, xend = 4.05, size = 0.25) +
  annotate(geom = "text", y = -0.45, x = 6.5, hjust = 0.5, vjust = 0, label = "Mean estimated change", size = 3) +
  annotate(geom = "segment", y = -0.5, yend = -0.3, x = 6.65, xend = 6.9, size = 0.25)

# % < 0

m1.fd_out0 <- apply(m1.fd, 2, function(x) ifelse(median(x) < 0, sum(x < 0) / length(x), sum(x > 0) / length(x)))
m1.fd_outROPE <- apply(m1.fd, 2, function(x) ifelse(median(x) < 0, sum(x < 0 - ROPE) / length(x), sum(x > 0 + ROPE) / length(x)))

fd_changelab <- ifelse(apply(m1.fd, 2, median) < 0, "Decrease by", "Increase by")
fd_changelab_short <- ifelse(apply(m1.fd, 2, median) < 0, "-", "+")
fd_changetext <- paste(fd_changelab, abs(round(apply(m1.fd, 2, median) * 100)), "points", sep = " ")
fd_changetext_short <- paste(fd_changelab_short, abs(round(apply(m1.fd, 2, median) * 100)), sep = "")

m1.fd_annotate <- data.frame(x = apply(m1.fd, 2, median), 
                             y = rev(unique(fd_long$Varpos)), 
                             out0 = paste(round(m1.fd_out0 * 100, digits = 1), "%", sep = ""), 
                             outROPE = paste(round(m1.fd_outROPE * 100, digits = 1), "%", sep = ""),
                             change = fd_changetext_short)
fd.p <- fd.p + geom_text(data = m1.fd_annotate, aes(y = x, x = y, label = outROPE), color = "white", nudge_x = 0.075, size = 3) + 
  geom_text(data = m1.fd_annotate, aes(y = x, x = y, label = change), color = "black", nudge_x = -0.1, size = 4)

ggsave(fd.p, file = "Output_Tables-and-Figures/crises_nomult_hiviol_m1_fd.pdf", width = 6, height = 8)